from src.SDRmodel import SDR_model
from src.ShallowDonor import Shallow_Donor
from src.Bulk import Silicon
import numpy as np

#import matplotlib.mlab as mmlab
import matplotlib.pyplot as plt

#--------------------------------------
#-------------- To Do list ------------
# 1) spin relaxation:
#	a- donor
#	b- electron
#	c- holes
#--------------------------------------

np.set_printoptions(suppress=True)

Bi = Shallow_Donor("Bi")
P = Shallow_Donor("P")
Si = Silicon()

#specie,d=0,g=0,c=0,rp=0,rap=1,w=0
SDR = SDR_model(Si, Bi, rp=0.01, rap=0.1, w=1, temperature=16, average_dist=1)

B = 0.1# Tu ne la mets pas en Gauss de merde

SDR.NRrates = SDR.Perform_NRrates(B)
NRpop, norm = SDR.Perform_population(SDR.NRrates, False)
# transitions =[[M1,b1,M2,b2],[M1,b1,M2,b2],...]
#NMR
#transitions=[[5,1,4,1],[4,1,3,1],[3,1,2,1],[2,1,1,1],[1,1,0,1],[0,1,-1,1],[-1,1,-2,1],[-2,1,-3,1],[-3,1,-4,1],[-4,1,-5,1],[4,-1,3,-1],[3,-1,2,-1],[2,-1,1,-1],[1,-1,0,-1],[0,-1,-1,-1],[-1,-1,-2,-1],[-2,-1,-3,-1],[-3,-1,-4,-1]]
#EPR
SpectroXband_X = [0.0659, 0.0795, 0.0989, 0.1263, 0.165, 0.2178, 0.2851, 0.364, 0.4509, 0.5434]
SpectroXband_Y = [1, 0.3771, 0.204, 0.126, 0.08, 0.056, 0.101, 0.086, 0.112, 0.241]
transitions = [[5, 1, 4, -1, 0.0659], [4, 1, 3, -1, 0.0795], [3, 1, 2, -1, 0.0989], [2, 1, 1, -1, 0.1263], [1, 1, 0, -1, 0.165], [0, 1, -1, -1, 0.2178], [-1, 1, -2, -1, 0.2851], [-2, 1, -3, -1, 0.364], [-3, 1, -4, -1, 0.4509], [-4, 1, -5, 1, 0.5434]]
#EPR with B
#transitions=[[5,1,4,-1,B],[4,1,3,-1,B],[3,1,2,-1,B],[2,1,1,-1,B],[1,1,0,-1,B],[0,1,-1,-1,B],[-1,1,-2,-1,B],[-2,1,-3,-1,B],[-3,1,-4,-1,B],[-4,1,-5,1,B]]
#EPR
#transitions=[[5,1,4,-1],[4,1,3,-1],[3,1,2,-1],[2,1,1,-1],[1,1,0,-1],[0,1,-1,-1],[-1,1,-2,-1],[-2,1,-3,-1],[-3,1,-4,-1],[-4,1,-5,1]]
T0 = [4, 1, 3, -1]
Simu_Y = []
Simu_X = []
for T in transitions:
	B = T[4]
	#SDR.NRrates=SDR.Perform_NRrates(B)
	#-----------------Construction of the rate matrix
	SDR.Rrates = SDR.Apply_Transition(SDR.NRrates, T, B)
	#SDR.Rrates=SDR.Apply_Transition(SDR.Rrates,T0,B)
	#-----------------Solve the system
	Rpop, trash = SDR.Perform_population(SDR.Rrates, True, norm)
	#-----------------Calculate the number of recombining pairs
	#rp = SDR.Perform_Recombining_Pairs(B)
	#-----------------Calculate the number of carriers
	NoC = SDR.Perform_number_of_carriers(NRpop, Rpop)
	#-----------------Store the results
	Simu_Y.append(NoC)
	Simu_X.append(B)
Simu_Y = Simu_Y / np.max(Simu_Y)

#-----------------Compute the variance of the error
error = [Simu_Y[i] - SpectroXband_Y[i] for i in range(np.size(Simu_Y))]
error_var = np.var(error)
print "Error variance= ", error_var

print Simu_Y

l = plt.plot(Simu_X, Simu_Y, 'ro', SpectroXband_X, SpectroXband_Y, "bo")

plt.xlabel('Magnetic Field')
plt.ylabel('SDR intensity')
plt.title(r'$\mathrm{Simulation\ of\ SDR\ intensity:}\ $')#\mu=100,\ \sigma=15 
margin = 0.1
xmin = np.min(Simu_X) - margin * (np.max(Simu_X) - np.min(Simu_X))
xmax = np.max(Simu_X) + margin * (np.max(Simu_X) - np.min(Simu_X))
ymin = np.min((Simu_Y, SpectroXband_Y)) - margin * (np.max(Simu_Y) - np.min(Simu_Y))
ymax = np.max((Simu_Y, SpectroXband_Y)) + margin * (np.max(Simu_Y) - np.min(Simu_Y))
plt.axis([xmin, xmax, ymin, ymax])
plt.grid(True)

plt.show()
